2. Seminar 1. EEG analysis

Plan

1. Read and preprocess the data

2. Use ICA for noise reduction

3. Compute ERP

4. Plot topomaps for ERP

5. Compute beta band envelopes for ERP

6. Compute coherence

# For Colab only
#!pip install mne
# !pip install mne-connectivity

# !wget "https://www.physionet.org/files/eegmmidb/1.0.0/S003/S003R03.edf"

# !wget "https://www.physionet.org/files/eegmmidb/1.0.0/S003/S003R03.edf.event"

# !ls

import mne
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns
import mne_connectivity

C:\Users\Alexander\anaconda3\lib\site-packages\pandas\core\computation\expressions.py:21: UserWarning: Pandas requires version '2.8.4' or newer of 'numexpr' (version '2.7.3' currently installed).
  from pandas.core.computation.check import NUMEXPR_INSTALLED
C:\Users\Alexander\anaconda3\lib\site-packages\pandas\core\arrays\masked.py:60: UserWarning: Pandas requires version '1.3.6' or newer of 'bottleneck' (version '1.3.2' currently installed).
  from pandas.core import (

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
C:\Users\ALEXAN~1\AppData\Local\Temp/ipykernel_4836/4132493247.py in <module>
      5 import matplotlib
      6 import seaborn as sns
----> 7 import mne_connectivity

ModuleNotFoundError: No module named 'mne_connectivity'

%matplotlib notebook

mne.io includes the funtions for different EEG-record formats

Optional To Do File formats (EDF, FIFF)

sample = mne.io.read_raw_edf('S003R03.edf', verbose=False, preload=True)

Get some info about a record

sample.info

Measurement date	August 12, 2009 16:15:00 GMT
Experimenter	Unknown
Participant	Unknown
Digitized points	Not available
Good channels	64 EEG
Bad channels	None
EOG channels	Not available
ECG channels	Not available
Sampling frequency	160.00 Hz
Highpass	0.00 Hz
Lowpass	80.00 Hz

# Sampling frequency
sample.info['sfreq']

160.0

# Length in seconds
len(sample) / sample.info['sfreq']

125.0

# Number of channels
len(sample.ch_names)

64

2.1. Channel selection and adding a montage

sample.ch_names[:3]

['Fc5.', 'Fc3.', 'Fc1.']

# fix trailing dots in channel names
# use sample.rename_channels(map)

# YOUR CODE HERE
sample.rename_channels(lambda x: x.replace('.', ''))

Measurement date	August 12, 2009 16:15:00 GMT
Experimenter	Unknown
Participant	Unknown
Digitized points	Not available
Good channels	64 EEG
Bad channels	None
EOG channels	Not available
ECG channels	Not available
Sampling frequency	160.00 Hz
Highpass	0.00 Hz
Lowpass	80.00 Hz
Filenames	S003R03.edf
Duration	00:02:05 (HH:MM:SS)

sample.ch_names[:3]

['Fc5', 'Fc3', 'Fc1']

sample.ch_names

['Fc5',
 'Fc3',
 'Fc1',
 'Fcz',
 'Fc2',
 'Fc4',
 'Fc6',
 'C5',
 'C3',
 'C1',
 'Cz',
 'C2',
 'C4',
 'C6',
 'Cp5',
 'Cp3',
 'Cp1',
 'Cpz',
 'Cp2',
 'Cp4',
 'Cp6',
 'Fp1',
 'Fpz',
 'Fp2',
 'Af7',
 'Af3',
 'Afz',
 'Af4',
 'Af8',
 'F7',
 'F5',
 'F3',
 'F1',
 'Fz',
 'F2',
 'F4',
 'F6',
 'F8',
 'Ft7',
 'Ft8',
 'T7',
 'T8',
 'T9',
 'T10',
 'Tp7',
 'Tp8',
 'P7',
 'P5',
 'P3',
 'P1',
 'Pz',
 'P2',
 'P4',
 'P6',
 'P8',
 'Po7',
 'Po3',
 'Poz',
 'Po4',
 'Po8',
 'O1',
 'Oz',
 'O2',
 'Iz']

# 19 channels from 10-20. no A1 and A2 here
# Be careful. Pure 10-20 labeling differs from high-resolution montages
# In MNE, 10-20 montage is actually an extended high-resulution version of 10-20
# FYI, mapping from pure 10-20 to high-resolution versions
# T3 = T7
# T4 = T8
# T5 = P7
# T6 = P8

channels_to_use = [
    # prefrontal
    'Fp1',
    'Fp2',
    # frontal
    'F7',
    'F3',
    'F4',
    'Fz',
    'F8',
    # central and temporal
    'T7',
    'C3',
    'Cz',
    'C4',
    'T8',
    # parietal
    'P7',
    'P3',
    'Pz',
    'P4',
    'P8',
    # occipital
    'O1',
    'O2',
]

sample_1020 = sample.copy().pick_channels(channels_to_use)

# check that everything is OK
assert len(channels_to_use) == len(sample_1020.ch_names)

ch_map = {ch.lower(): ch for ch in sample_1020.ch_names}

ten_twenty_montage = mne.channels.make_standard_montage('standard_1020')
len(ten_twenty_montage.ch_names)

94

ten_twenty_montage.ch_names = [ch_map[ch.lower()] if ch.lower() in ch_map else ch 
                               for ch in ten_twenty_montage.ch_names]

sample_1020.set_montage(ten_twenty_montage)

Measurement date	August 12, 2009 16:15:00 GMT
Experimenter	Unknown
Participant	Unknown
Digitized points	22 points
Good channels	19 EEG
Bad channels	None
EOG channels	Not available
ECG channels	Not available
Sampling frequency	160.00 Hz
Highpass	0.00 Hz
Lowpass	80.00 Hz
Filenames	S003R03.edf
Duration	00:02:05 (HH:MM:SS)

sample_1020.plot_sensors(show_names=True);

2.2. Explore the signals

sample_1020.plot_psd();

NOTE: plot_psd() is a legacy function. New code should use .compute_psd().plot().
Effective window size : 12.800 (s)

2.2.1. Band-pass filtering

It’s better to remove low-freq components < 1 Hz and high-freq > 50Hz (non-informative for EEG)

Let’s use 4-th order Butterworth filter (default IIR filter)

sample_1020.filter(l_freq=1, h_freq=50, method='iir')

Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 1 - 50 Hz

IIR filter parameters
---------------------
Butterworth bandpass zero-phase (two-pass forward and reverse) non-causal filter:
- Filter order 16 (effective, after forward-backward)
- Cutoffs at 1.00, 50.00 Hz: -6.02, -6.02 dB

Measurement date	August 12, 2009 16:15:00 GMT
Experimenter	Unknown
Participant	Unknown
Digitized points	22 points
Good channels	19 EEG
Bad channels	None
EOG channels	Not available
ECG channels	Not available
Sampling frequency	160.00 Hz
Highpass	1.00 Hz
Lowpass	50.00 Hz
Filenames	S003R03.edf
Duration	00:02:05 (HH:MM:SS)

# Plot psd after filtering

# YOUR CODE HERE

2.2.2. Plot EEG signals

sample_1020.plot(n_channels=8, duration=20);

Using matplotlib as 2D backend.

# Plot in better scale. Use 'scaling' argument

# YOUR CODE HERE

2.3. Extracting events

Mne has several functions for event selection.

• mne.find_events is used when events are stored in trigger channels (e.g. FIFF format)

• mne.events_from_annotations is used for when events are stored in annotations (EDF+ format)

Look for documentation for your EEG-record format

Here we have EDF+ format

events, events_dict = mne.events_from_annotations(sample_1020)

Used Annotations descriptions: ['T0', 'T1', 'T2']

events_dict

{'T0': 1, 'T1': 2, 'T2': 3}

events[:5]

array([[   0,    0,    1],
       [ 672,    0,    3],
       [1328,    0,    1],
       [2000,    0,    2],
       [2656,    0,    1]])

epochs = mne.Epochs(sample_1020, events,  tmin=-0.5, tmax=0.8, preload=True)

Not setting metadata
30 matching events found
Setting baseline interval to [-0.5, 0.0] sec
Applying baseline correction (mode: mean)
0 projection items activated
Using data from preloaded Raw for 30 events and 209 original time points ...
1 bad epochs dropped

pd.DataFrame(epochs.events, columns=['_', '__', 'event_id'])['event_id'].value_counts()

1    14
3     8
2     7
Name: event_id, dtype: int64

Check that length is right

for epoch in epochs:
    break
epoch.shape

(19, 209)

epoch.shape[1] / sample_1020.info['sfreq']

1.30625

sample_1020.to_data_frame().shape

(20000, 20)

df = epochs.to_data_frame()
df.head(3).iloc[:, :10]

	time	condition	epoch	C3	Cz	C4	Fp1	Fp2	F7	F3
0	-0.50000	3	1	-4.547919	-1.344287	1.476721	236.675891	244.275703	125.936913	88.829681
1	-0.49375	3	1	-4.844281	-12.002902	3.515891	175.970007	183.425016	103.803572	63.009365
2	-0.48750	3	1	-2.909761	-11.900140	6.341468	127.095181	132.985975	61.669835	31.478465

df[sample_1020.ch_names + ['epoch']].groupby('epoch').agg(lambda arr: arr.max() - arr.min()).hist(figsize=[10, 10]);
plt.tight_layout()

epochs = mne.Epochs(sample_1020, events,  tmin=-0.5, tmax=0.8, reject={'eeg': 600e-6}, preload=True, baseline=(-.1, 0))

Not setting metadata
30 matching events found
Applying baseline correction (mode: mean)
0 projection items activated
Using data from preloaded Raw for 30 events and 209 original time points ...
    Rejecting  epoch based on EEG : ['Fp1', 'Fp2']
    Rejecting  epoch based on EEG : ['Fp2']
    Rejecting  epoch based on EEG : ['Fp2']
    Rejecting  epoch based on EEG : ['Fp2']
    Rejecting  epoch based on EEG : ['Fp2']
    Rejecting  epoch based on EEG : ['Fp1', 'Fp2']
    Rejecting  epoch based on EEG : ['Fp1', 'Fp2']
    Rejecting  epoch based on EEG : ['Fp1', 'Fp2']
    Rejecting  epoch based on EEG : ['Fp1', 'Fp2']
    Rejecting  epoch based on EEG : ['Fp1', 'Fp2']
    Rejecting  epoch based on EEG : ['Fp2']
    Rejecting  epoch based on EEG : ['Fp1', 'Fp2']
    Rejecting  epoch based on EEG : ['Fp1', 'Fp2']
    Rejecting  epoch based on EEG : ['Fp1', 'Fp2']
15 bad epochs dropped

PSD on epochs differs from the raw. More averaging is used

epochs.plot_psd();

NOTE: plot_psd() is a legacy function. New code should use .compute_psd().plot().
    Using multitaper spectrum estimation with 7 DPSS windows
Averaging across epochs...

epochs.plot(n_channels=8, scalings={'eeg':1e-4});

epochs.event_id

{'1': 1, '2': 2, '3': 3}

# check number of events of each type
# use epochs.events

# Your code here

pd.DataFrame(epochs.events, columns=['_', '__', 'event_id'])['event_id'].value_counts()

3    6
1    5
2    4
Name: event_id, dtype: int64

evoked_T0 = epochs['1'].average()
evoked_T1 = epochs['2'].average()
evoked_T2 = epochs['3'].average()

evoked_T0.plot(spatial_colors=True);

evoked_T1.plot(spatial_colors=True);

evoked_T2.plot(spatial_colors=True);

2.4. Independent Component Analysis for Artifact Removal

ica = mne.preprocessing.ICA(n_components=10, random_state=42)

ica.fit(sample_1020)

Fitting ICA to data using 19 channels (please be patient, this may take a while)
Selecting by number: 10 components
Fitting ICA took 1.4s.

Method	fastica
Fit	17 iterations on raw data (20000 samples)
ICA components	10
Available PCA components	19
Channel types	eeg
ICA components marked for exclusion	—

ica.plot_sources(sample_1020);

Creating RawArray with float64 data, n_channels=10, n_times=20000
    Range : 0 ... 19999 =      0.000 ...   124.994 secs
Ready.

ica.plot_components();

Inspect ICA components more deeply. Check out spectrogram. Segments info is not very relevant here since we build ICA on the raw data

We expect to see alpha and beta rythms picks on the spectrogram for good components (7-13 Hz and 13-30Hz respectively). And also slight decrease as frequency goes higher

ica.plot_properties(sample_1020, picks=[4]);

    Using multitaper spectrum estimation with 7 DPSS windows
Not setting metadata
62 matching events found
No baseline correction applied
0 projection items activated

ica.plot_overlay(sample_1020, exclude=[0, 1, 4, 5, 8, 9], picks=['F3']);

Applying ICA to Raw instance
    Transforming to ICA space (10 components)
    Zeroing out 6 ICA components
    Projecting back using 19 PCA components

ica.exclude = [0, 1]

sample_1020_clr = sample_1020.copy()

ica.apply(sample_1020_clr)

Applying ICA to Raw instance
    Transforming to ICA space (10 components)
    Zeroing out 2 ICA components
    Projecting back using 19 PCA components

Measurement date	August 12, 2009 16:15:00 GMT
Experimenter	Unknown
Participant	Unknown
Digitized points	22 points
Good channels	19 EEG
Bad channels	None
EOG channels	Not available
ECG channels	Not available
Sampling frequency	160.00 Hz
Highpass	1.00 Hz
Lowpass	50.00 Hz
Filenames	S003R03.edf
Duration	00:02:05 (HH:MM:SS)

# plot channels

# YOUR CODE HERE

epochs = mne.Epochs(sample_1020_clr, events,  tmin=-0.5, tmax=0.8, reject={'eeg': 600e-6}, preload=True, baseline=(-.1, 0))

Not setting metadata
30 matching events found
Applying baseline correction (mode: mean)
0 projection items activated
Using data from preloaded Raw for 30 events and 209 original time points ...
1 bad epochs dropped

evoked_T0 = epochs['1'].average()
evoked_T1 = epochs['2'].average()
evoked_T2 = epochs['3'].average()

evoked_T0.plot(spatial_colors=True);

evoked_T1.plot(spatial_colors=True);

evoked_T2.plot(spatial_colors=True);

evoked_T0.plot_topomap(times=[0, .2, .4, .6, .8], vlim=(-50,50));

evoked_T1.plot_topomap(times=[0, .2, .4, .6, .8], vlim=(-50,50));

evoked_T2.plot_topomap(times=[0, .2, .4, .6, .8], vlim=(-50,50));

2.5. Dynamics of alpha and beta activity

evoked_T0_alpha = evoked_T0.copy().filter(l_freq=7, h_freq=13, method='iir', verbose=False).apply_hilbert(envelope=True)
evoked_T1_alpha = evoked_T1.copy().filter(l_freq=7, h_freq=13, method='iir', verbose=False).apply_hilbert(envelope=True)
evoked_T2_alpha = evoked_T2.copy().filter(l_freq=7, h_freq=13, method='iir', verbose=False).apply_hilbert(envelope=True)

evoked_T1_alpha.plot(spatial_colors=True);

evoked_T0_alpha.plot_topomap(times=[0, .1, .2, .3, .4, .6], vlim=(0,30));

evoked_T1_alpha.plot_topomap(times=[0, .1, .2, .3, .4, .6], vlim=(0,30));

evoked_T2_alpha.plot_topomap(times=[0, .1, .2, .3, .4, .6], vlim=(0,30));

evoked_T0_beta_low = evoked_T0.copy().filter(l_freq=13, h_freq=20, method='iir', verbose=False).apply_hilbert(envelope=True)
evoked_T1_beta_low = evoked_T1.copy().filter(l_freq=13, h_freq=20, method='iir', verbose=False).apply_hilbert(envelope=True)
evoked_T2_beta_low = evoked_T2.copy().filter(l_freq=13, h_freq=20, method='iir', verbose=False).apply_hilbert(envelope=True)

evoked_T0_beta_low.plot_topomap(times=[0, .2, .4, .6, .8], vlim=(0,30));

evoked_T1_beta_low.plot_topomap(times=[0, .2, .4, .6, .8], vlim=(0,30));

evoked_T2_beta_low.plot_topomap(times=[0, .2, .4, .6, .8], vlim=(0,30));

2.6. Computing functional connectivity

conn_T1, freqs, times, _, _ = mne_connectivity.spectral_connectivity(epochs['2'], method='coh')

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
C:\Users\ALEXAN~1\AppData\Local\Temp/ipykernel_7908/3959819980.py in <module>
----> 1 import mne_connectivity
      2 conn_T1, freqs, times, _, _ = mne_connectivity.spectral_connectivity(epochs['2'], method='coh')

ModuleNotFoundError: No module named 'mne_connectivity'

def plot_topomap_connectivity_2d(info, con, picks=None, pairs=None, vmin=None, vmax=None, cm=None, show_values=False, show_names=True):
    """
    Plots connectivity-like data in 2d
    
    Drawing every pair of channels will likely make a mess
    There are two options to avoid it:
    - provide picks
    - provide specific pairs of channels to draw
    """
    
    # get positions
    _, pos, _, _, _, _, _ = mne.viz.topomap._prepare_topomap_plot(info, 'eeg');
    
#     if picks is None and pairs is None:
#         picks = info.ch_names
    
    ch_names_lower = [ch.lower() for ch in info.ch_names]
    if picks:
        picks_lower = [ch.lower() for ch in picks]
    if pairs:
        pairs_lower = [tuple(sorted([ch1.lower(), ch2.lower()])) for ch1, ch2 in pairs]
    
    rows = []
    for idx1, ch1 in enumerate(ch_names_lower):
        for idx2, ch2 in enumerate(ch_names_lower):
            if ch1 >= ch2:
                continue
            if picks and (ch1 not in picks_lower or ch2 not in picks_lower):
                    continue
            if pairs and (ch1, ch2) not in pairs_lower:
                    continue
            rows.append((
                pos[idx1],
                pos[idx2],
                con[idx1, idx2]
            ))
    
    if not len(rows):
        raise ValueError('No pairs to plot')
    
    con_to_plot = np.array([row[2] for row in rows])
    if vmin is None:
        vmin = np.percentile(con_to_plot, 2)
    if vmax is None:
        vmax = np.percentile(con_to_plot, 98)
    norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
    
    if cm is None:
        cm = sns.diverging_palette(240, 10, as_cmap=True)
    
    fig, ax = plt.subplots(figsize=[5, 5])
    mne.viz.utils.plot_sensors(info, show_names=show_names, show=False, axes=ax);
    for row in rows:
        rgba_color = cm(norm(row[2]))
        plt.plot([row[0][0], row[1][0]], [row[0][1], row[1][1]], color=rgba_color)
        if show_values:
            plt.text((row[0][0] + row[1][0]) / 2, 
                     (row[0][1] + row[1][1]) / 2, 
                     '{:.2f}'.format(row[2]))

conn_T0, freqs, times, _, _ = mne_connectivity.spectral_connectivity(epochs['1'], method='coh', verbose=False);
conn_T1, freqs, times, _, _ = mne_connectivity.spectral_connectivity(epochs['2'], method='coh', verbose=False);
conn_T2, freqs, times, _, _ = mne_connectivity.spectral_connectivity(epochs['3'], method='coh', verbose=False);

freqs

conn_T0_beta = conn_T0[:, :, 12:27].mean(axis=2)
conn_T0_beta = conn_T0_beta + conn_T0_beta.T

conn_T1_beta = conn_T1[:, :, 12:27].mean(axis=2)
conn_T1_beta = conn_T1_beta + conn_T1_beta.T

conn_T2_beta = conn_T2[:, :, 12:27].mean(axis=2)
conn_T2_beta = conn_T2_beta + conn_T2_beta.T

plot_topomap_connectivity_2d(epochs.info, conn_T1_beta, picks=epochs.ch_names);

plot_topomap_connectivity_2d(epochs.info, conn_T0_beta, 
                             pairs=[('F7', 'F4'), ('O2', 'T7'), ('C3', 'C4'), ('P7', 'P8'), ('F8', 'T8'), ('O1', 'O2'), ('O1', 'P4')],
                             show_values=True,
                             show_names=False
                            
                            );

# calculate coherence in alpha band
# plot 5-10 pairs that you are interested in

# YOUR CODE HERE